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ABSTRACT 

A numerical simulation model for 1-D steady state flow in open channels with ir- 
regular non-uniform cross section and a numerical simulation model for 2-D pollu- 
tant transport using a modified form of the mixing analysis based on the concept 
of stream tubes(MABOCOST) are presented. Gradually varied flow(GVF) in the 
channel is considered. The flow simulation model is formulated using the fourth order 
Runge Kutta method. The flow model is solved for an open channel with specified 
irregular cross sections assuming immobile bed conditions.The transport simulation 
model utilizes the flow simulation solutions to simulate pollutant transport in the 
channel, assuming advective and dispersive transport processes.The transport model 
simulates the spatial and temporal distribution of pollutant concentration for given 
pollution sources. The simulated results are compared with reported experimental 
results for a similar 2-D system as well as reported 2-D simulation results using 
MABOCOST. These comparisions and mass balance analysis show that the com- 
bined 1-D flow and 2-D pollutant transport simulation results may be used as ap- 
proximate solutions for the 2-D system for preliminary assessment. Both conser^atiie 
pollutants and pollutants undergoing decay process are considered. 
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cross-sectional area of discharge (m^) 
depth- averaged concentration {mg/L) 
dispersion factor (s~^) 
longitudinal mixing coefficient {m^/s) 
transverse mixing coefficient (m^/s) 
accelaration due to gravity (m/s^)) 
flow depth of stream channel (m) 

ratio between distance travelled in time step to longitudinal length of grid element 

fraction of overlapping area of element with its nth adjacent element on left- 

longitudinal and transverse metric coefficients respectively in stream-tube equation 

Manning’s coefficient of roughness 

normalized cumulative discharge (^) 

flow rate of stream channel {w? js) 

cumulative discharge (m®/s) 

hydraulic radius of the stream channel (m) 

fraction of overlapping area of element with its nth adjacent element on right 



bottom slope of the stream channel 
energy slope of the channel flow 
top width of flow cross-section (m) 
time (s) 

depth-averaged longitudinal velocity (m/s) 

longitudinal coordinate axis (m) 

depth of flow with respect to a datum (m) 

transverse coordinate axis (m) 

rate of growth or decay of pollutant (s“^) 

magnitudeof source-sink term (mp/L)(L/s) 

concentration surplus between element and its neighbor, (mg/L) 

stream-tube element volume m^) 



\'m 


LIST OF FIGURES 


Figure Title Page 

2.1 Schematic Diagram of the Stream Tubes 10 

2.2 Division of Stream Tubes into variable length elements 10 

3.1 Channel Sections IS 

3.2 Comparison of Simulation Results with Experimental and 22 

original MABOCOST Results 

3.3 Concentration profiles for all Stream Tubes at X = Am 26 

3.4 Comparison of Simulation Concentrations at X = 8m 29 

for different a 

3.5 Concentration profiles for multiple Slug Injections 30 



IX 


LIST OF TABLES 


Table 

Title 

Page 

3.1 

Flow parameter values 

17 

3.2 

Pollutant transport parameter values 

19 

3.3 

Mass balance calculations 

25 


Chapter 1 


Introduction 


1.1 Introduction 

The human interference in the natural activities can be seen through so many envi- 
ronmental hazards. The manifestation of such an interference is poor water quality 
in natural Rivers which is ultimately causing great problems in maintaining ecolog- 
ical balance. As such it has long been noted that the problem of poor water quality 
in the natural rivers interferes with the human life and the manifestation of that 
interference is indicated by disease transmission and aesthetic nuisances. Nowadays, 
water quality in natural sreams is becoming a major issue all over the world due to 
wide spread pollution. Hence, the necessity of evolving effective remedial measures 
are being very important. 

On the contrary, modelling the transport of pollutants is difficult due to nature's 
hetrogeneity and the complexity of estimating the flow and transport parameters 
especially, in natural flow systems like meanders, irregular cross sections, tribu- 
taries and lateral flows. Due to these complexities, an effective control of pollution 
and maintaining aesthetic conditions in natural streams requires development and 
solution of water quality models taking into consideration the complete processes 
involved. 

In this study (i) A numerical model is developed for simulating 1-D flow- under 
steady state conditions for channels with irregular cross sections (ii) A simulation 
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model is solved for simulating the 2-D transport of a pollutant in the channel. This 
model is applicable for both 1-D and 2-D flow systems using the stream tube concept. 

The transport model simulates the spatial and temporal distribution of pollutant 
concentration for given pollution sources. The simulatoion results for 1-D simula- 
tion are compared with reported experimental results for a similar 2-D system as 
well as reported 2-D simulation results using MABOCOST. These comparisions and 
mass balance analysis show that the 1-D simulation results can be used as approx- 
imate solutions for the 2-D system for preliminary assessment. Both conservative 
pollutants and pollutants undergoing decay process are considered. 

1.2 Literature Review 

A steady non-uniform flow in a channel with gradual changes in its water surface 
ele\’’ation is termed as gradually- varied flow (GVF). In GVF, the depth of flow and 
velocity vary along the channel. The bed slope, water surface slope, and energy slope 
will all differ from each other. Because of its practical importance the computation 
of GVF has been a topic of continued interest for somany years. 

To meet the practical needs, various solution procedures involving graphical and 
numerical methods were evolved. The advent of high speed computers has given rise 
to general programmes utilizing sophisticated numerical techniques for solving GVF 
in natural channels. The various available procedures for computing GVF profiles 
can be classified as: 

1. Direct integration 

2. Numerical method 

3. Graphical method 

Computation procedures applicable to an artificial channel may not be directly 
applicable to natural channels (Subramanya,1989) of irregular cross section since in 
a natural channel, the cross sectional properties are known only at specified locations 
while in prismatic artificial channel the cross-sections are generally constant all along 
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the channel. Because of this, some methods have been developed, particularly for 
use in natural channels, e.g. Standard-step metod. 

The numerical solution procedures to solve GVF problems can be broadly clas- 
sified into two categories as: 

(I) Simple Numerical Methods 

(II) Advanced Numerical Methods 

Simple numerical methods usually attempt to solve the energy equation either in 
the form of the differential energy equation of GVF or in the form of the Bernouli 
equation. Two commonly used simple numerical methods to solve GVF problems 
are as follows: 

(i) Direct-step method and 

(ii) Standard step method 

Advanced numerical methods are normally suitable for use in digital computers as 
they involve a large number of repeated calculations. The basic differential equation 
of GVF can be expressed as a function of depth of flow only for a given bed slope, 
Manning’s coefficient of roughness, flow and channel geometry. A class of methods 
which is particularly suitable for numerical solution of this non-linear equation is 
the Runge-Kutta method. There are different orders of Runge-Kutta methods and 
all of them evaluate the depth of flow at a desired section provided the depth at a 
known section is given. The fourth order method is used here for its accuracy. The 
various numerical methods are as follows: 

(a) Standard Fourth Order Runge-Kutta Method (SRK) 

(b) Kutta-Merson Method (KM) 

(c) Trapezoidal Method (TRAP) 

Studies have, been reported (in Subramanya,1989) by Apelt (1971) and Humpridge 
and Moss (1972) on the SRK method;by Apelt on the KM and TRAP methods and 
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by Prasad (1970) on the TRAP method. It has been found that all these three 
methods are capable of direct determination of the GVF profiles both in upstream 
and down stream directions irrespective of the nature of flow, i.e whether the flow 
is subcritical or super critical. It has been observed that the SRK and KM methods 
require less computational effort than TRAP method. Also, while the SRK method 
is slightly more efficient than the KM method, the possibility of providing automatic 
control of the step size and accuracy in the KM process makes it a strong contender. 
Though there are host of Graphical methods, these have become obsolete with the 
advent of easy access to digital computers. In this study SRK method is adopted 
for computation. 

The Preissmann scheme (Chaudhary,1990) has been extensively used for the 
past three decades. It has the advantages that a variable spatial grid may be used. 
By varying the weighting coefficient steep wave fronts may be properly simulated. 
Preissmann scheme yields an exact solution of the linearized form of the governing 
equation for a particular value of step size in both the directions of length and time. 

For N number of reaches on the channel this scheme yields 2N equations for each 
grid point. The equations cannot be written for the grid points at the downstream 
end. However, the number of unknowns are 2{N + 1), that is, two unknowns for 
each grid point. Thus, two more equations are needed for a unique solution. These 
are provided by the end conditions. Unlike the explicit schemes, the equations 
describing the end conditions are directly included in the system of equations. In 
other words, we do not have to use the characteristic equations or the reflection 
procedures. Because steady state conditions are assumed in this solution, the SRK 
method is adopted. 

The scheme is unconditionally stable provided the flow variables are weighted 
towards the next higher time level. An unconditional stability means that there is 
no restriction on the step size of both the length and time. However, for accuracy 
the courant number should be close to 1.0 for frictionless channel. But if linearized 
form , of the friction loss term is included in the analysis then Vedernikov number 
should be close to 1.0 
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As far as pollution transport in streams is concerned, mathematical models of 
different le\els of sophistication, based on some finite difference representations of 
mass balance equation, ha\’'e been developed. The mixing model is based either on 
the rectangular channel assumption or the stream tube concept (Beitaos,1979; Paily 
and Sayre,1978; Rood and Holley, 1974; Stefan and Gulliver, 1978; Yotsukura,1977; 
Yotsukura and Cobb,1972). The latter class of models have become more popular 
in recent years because of their ability to simulate processes in natural stream chan- 
nels having irregular cross sections. Numerical models are far more advantageous 
than analytical models, with respect to flexibility, as they provide the capability 
of accounting for actual flow patterns and boundary geometry as well as spatially 
varying mixing parameters. 

Most of the mixing models developed to-date are one-dimensional in nature, that 
is, complete mixing of pollutants across the stream sections is assumed. Models 
given by Fukuoka and Sayre (1973), Valentine and Wood (1977), Beltaos(1980), 
Terragni and Young (1983) and McBribe and Rutherford (1984) are such examples. 
However, Elhadi (1984) in his report stated that except for very small streams, it 
has been shown that complete sectional mixing is not achieved until the pollutant 
has travelled long distances downstream from the point of release. 

The Two-dimensional mixing models developed to-date can be classified accord- 
ing to whether pollutant inputs are steady or unsteady. For example, McCorquodale 
et al.(1983) developed a mixing model for straight channels and models of Lau and 
Krishnappan (1981), Somlyody(1982) and Gowda(1984) utilized the stream-tube 
concept to represent non-uniform, meandering shapes of Rivers. Unsteady-state 
mixing models include those developed by Verboom (1974), Holly (1975) and On- 
ishi(1981) which deal with straight channels. In addition Holly and Gunge (1975) 
and Luk et al (1990) developed an unsteady-state stream-tube dispersion model. 
Also, Harden and Shen (1979) developed a transient mixing model based on a com- 
bined implicit, explicit finite difference scheme. 
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1.3 Objective 

This study describes a method qf modeling the two-dimensional mixing of pollutants 
in natural channels, with the use of a simple verification using analytical solutions 
and a reputed numerical method. The following are the objectives of this study. 

1. A 1-D flow simulation model for steady state GVF in an open channel is de- 
veloped and numerically solved using the Fourth order Runge-Kutta method. 

2. A modified form of the MABOCOST model is solved for 1-D flow conditions, 
with identical lengths of stream tubes at a given cross section, while considering 
2-D transverse mixing effects through the transverse mixing coefficient. 

3. The pollutant transport is simulated for a straight channel length injected 
pollutant slugs for conservative and decaying pollutants. 

4. The simulated pollutant concentrations for the idealized straight channel con- 
figuration is compared with reported experimental and numerical results for 
2-D flow in a mildly curved channel with similar cross sections as assumed in 
this study. 

5. The simultion results are also validated using the mass balance criterion for 
the pollutants. 



Chapter 2 


Methodology 


2.1 Introduction 

The flow simulation model considers one-dimensional Gradually Varied Flow and 
irregular cross sections of the channel. The equation of G\T is non-linear in nature. 
The Standard Fourth Order Runge-Kutta Method (SRK) is used for the solution of 
the flow equations for steady state conditions and immobile boundary. 

The pollutant transport equation is solved using a modified version of the Mixing 
Analysis Based On the Concept Of Stream Tubes (MABOCOST) using the time- 
fractional step method (Luk et al,1990). The governing equations for flow and 
pollutant transport along with the method of solution used are presented in this 
chapter. 

2.2 Governing Equations for Open Channel Flow 

The basic differential equation of GVF is derived from the concept of total energ}- 
in open channels. Since the water surface, in general, varies in the longitudinal 
(x) direction, the depth of flow and total energy are functions of x. The two basic 
assumptions in the derivation of Gradually Varied Flow equation are: 

1. The pressure distribution at any section is assumed to be hydrostatic. 
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2. At any depth the resistance to flow is assumed to be given by corresponding 
uniform-flow equation. If Manning’s formula is used, the slope term to be used 
is energy slope but not the bed slope. 


The differential equation of GVF also known as the dynamic equation of GVF 
(Subramanya,1989) is as follows: 


dx 



1 

^ S.43 


( 2 . 1 ) 


where, 

y = depth of flow with respect to a datum 
X = longitudinal coordinate 

5o = the bottom slope of the channel 
Sf = the energj" slope of the channel flow 
a = the kinetic-energy correction factor 

Q = the discharge in the channel 

T = the top width of flow cross section 
g = accelaration due to gravity 

A = cross sectional area of discharge 

The energy slope is defined as : 


( 2 . 2 ) 


2.3 Governing Equations for Pollutant Transport 
in Open Channels 

In the pollutant transport model the convection can be deflned as the transfer of a 
pollutant from one element in a stream tube to the corresponding element down- 
stream without any change in the total mass. The process of distribution of the 
pollutant to the adjacent elements without any change in the total mass of the 
pollutant concentration is known as dispersion. The decay causes reduction in the 
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total mass of the concentration in each element. In the transport model the pollu- 
tant dispersion is governed by the principle of conservation of mass. The continuity 
equation for neutrally buoyant substance can be written using natural coordinate 
system (Yotsukura and Sayre, 1976) as : 


dt rrix dx dn V <3^ 


-b aC -f P 


(2.3) 


This equation the case where the pollutant undergoes a first order chemical reaction, 
and where either a source or a sink is present owing to localized chemical or biological 
reactions. 


where, 

C = depth-averaged pollutant concentration 

u = depth-averaged longitudinal velocity 

h = local flow depth 

Ez = turbulent lateral transverse mixing coefficient 
Q = flow rate 

a = rate of chemical reactiion 

P = magnitude of the source or sink 

X = longitudinal coordinate 

q = cumulative discharge 

n = The non dimensional cumulative discharge 

rria: and niz are the metric coefficients to account for the length variation between 
coordinate surfaces (for curved channels nix and rriz are greater than 1.0 , for straight 
channels = mz=1.0). 

The stream tube coordinate system and the schamatic diagram showing stream 
tubes divided into varible length elements are shown in Fig 2.1 and Fig 2.2 respec- 
tively (Luk et al,1990). The nondimensional cumulative discharge ’n’ is the fraction 
of total flow between the bank and the lateral distance z ,that is, 


Q 


Q 



(2.4) 
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The performance of numerical schemes depends on the courant number. In the 
numerical grid formulated for the finite difierence solution of the 2-D mass balance 
equation the Courant number, L, is defined as : 


L = 


uAt 

Ax 


(2.5) 


where, 

At = time step length and 
Ax = longitudinal grid size 


Djordjevic and Radojkovic (1986) independently concluded that one of the simplest 
numerical algorithms that ensures stability and accuracy is the first order, asymmet- 
rical, explicit scheme, with the courant number set to unity throughout the entire 
grid. 

The MABOCOST method is based on the concept of stream tubes. Beltaos 
(1978) suggested that each stream tube must be divided into the elements with the 
length of each individual element equal to the product of the local velocity and time 
interval to satisfy the condition of the Courant number set to unity. The length of 
each idividual element in the stream tube is defined as : 


Ax = uAt (2.6) 

In the Eqn.2.3, the longitudinal mixing coefficient is assumed to be negligible as in 
natural channels the transverse mixing coefficient is few orders of magnitude higher 
than the longitudinal mixing coefficient. 

2.4 Solution Methodology 

The solution procedures adopted in this study for solving the flow and pollutant 
transport equations are discussed in this section. The following assumptions are 
made for solving the flow equations : 

1. The flow in the stream is steady state. 
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2. The flow is one dimensional. 

3. The velocity is constant across a section at. any given location. 

2.4.1 Methodology for Solution of Flow Model 

The 1-D numerical model for open channel flow involves solving the flow equation 
(Eqn.2.1) using the fourth order Runge-Kutta method. 

The fourth order Runge-Kutta method for solving the non-linear differential 
equation of GVF in open channels is as follows. 

Vi+l = yi+ 0(-^l 2-^2 + SiFs + K^) (2.7) 

where, 



F'iy) = 

% 

dx 

(2.8) 


Ki = 

^xF{yi) 

(2.9) 


= 

^F(yi + ^) 

(2.10) 


Kz = 


(2.11) 


K, = 

AxF{yi + Kz) 

(2.12) 

and 




Vi 

= The depth of flow at cross 

section 


Ax - 

= The longitudinal step size 




Natural channels are, in general, irregular in cross section. In irregulr cross sectional 
channels the area and top width are functions of the depth of flow. For each depth 
of flow the area and top width can be obtained graphically or be computed. A table 
is prepared constituting the area, top width and the corresponding depth for input 
to the computer. This procedure has been followed at each known cross section for 
this study. 

The channel is divided into N parts of known length interval Ax. Initially, the 
area of flow and top width are interpolated from the input table corresponding to 
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given depth of flow. The value of Ki is evaluated using Eqn. 2.8 for the known 
depth of flow. Knowing Ki, K 2 ,K 3 and K 4 are computed using Eqns. 2.10 to 2.12 . 
Using these computed values of Kx^K^.Kz and K 4 and yi the depth j/j+i at the next 
cross section is obtained. This process is continued for subsequent cross sections. 
The top width and cross sectional area of flow computing to a given depth yi is 
again obtained by interpolation from the table of A and T versus y. The average 
velocity of flow is determined by using the cross sectional area of flow and discharge 
values at each cross section. 

2.4.2 Methodology for Solution of Pollutant Transport Model 

The two-dimensional numerical model for transient pollutant transport involves solv- 
ing the three computational modules representing advection, dispersion and decay 
in sequence. 

Verboom(1975) showed that the two dimensional mass balance equation (Eqn. 

2.3) can be solved by using time fractional-step method. The basic priniciple of 
this method is the simplification of the solution of a multidimensional mass balance 
equation by splitting it into many smaller and simpler steps. Stability and consis- 
tency of the entire solution can be assured provided the numerical scheme chosen 
for each indiviual step is stable and consistent. 

Some assumptions are made in this transport simulation model. They are : 

1. The length of each individual element in each stream tubes is the product of 
local velocity and time step (maintaining the Courant number as 1.0) 

2. The velocity u at any particular cross section as obtained from the solution of 
the 1-D flow model is uniform across the width. However, for pollutant trans- 
port purpose the transverse mixing coefficient is assumed to be non zero 
to account for actual transverse velocity variation and transverse dispersion. 

3. The numerical dispersion and diffusion is avoided with the construction of 
the grid in a way that the length of each element equal to the distance that 
pollutant travels in one time step, by advection. 
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4. The length of the elements in all stream tubes at any particular cross section 
are equal due to 1-D flow simulation. 

A steady reach is divided into a number of stream tubes. Each stream tube is 
divided into variable length elements at different cross sections. Thus, a numerical 
grid is formulated with number of stream tubes and elements. The convective por- 
tion of the mass balance equation is first solved. The solution is obtained simply by 
advecting the pollutant forward to the next element, that is 


C{iJ,t + At) = C{i,j-l,t) (2.13) 

where, 

i =z stream tube and 
j = element in the stream tube. 


Following this step, the newly convected element concentrations are dispersed lat- 
erally to each of the adjacent elements. This is represented by the following finite 
difference representation: 


^/. . V X • X At ^ (lAxihAzDAC), 

c(i,],t + At) = + — ii 


p=zl 


{lAxi)p\{hAzD)p - hAzD\ACp ^ {rAxrhAzDAC\ 

p=l 


An^ 


-f 


(rA 2 :,)p|(/iAzD)p - hAzD\ACp 
An2 


(2.14) 


Where, 

D = the dispersion factor given by, 


Q2 

9 = the volume of an element 

Ax = the element length 
Az = the element width 
AC = the concentration excess and 


(2.15) 
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An = the fraction of the cummulative discharge between the centers 
of two elements 

Ez = The transverse mixing coefficient 

l,r are the fractions of overlapping area between an element and its adjacent element 
on the left and right sides, respectively. The subscript p represents the quantities 
for the p^^ adjacent elements, the subscripts I and r represent the quantity for the 
left and right adjacent elements respectively. Overbars represent quantities averaged 
between an element and its adjacent elements. 

Next , the concentrations thus obtained are changed according to a first order 
chemical reaction or decay. 

C{i,j,t + At) = C{i,j,t).exp{aAt) (2.16) 

Finally the source- sink terms for the current time step can be added to each in- 
dividual element. These processes are repeated for all time steps to obtain the 
concentrations for all elements at every time step. 

The depth of flow and velocity for each element in a stream tube is calculated 
by linear interpolation of its depth and velocity obtained by using the SRK method 
at two sections between which the element is situated. Using these depth and 
velocities for each idividual element as the input for a given value of transverse 
mixing coefficient and a time step the dispersive part of the transport model is 
solved. 

The pollutant released in the elements at a given location is convected to the next 
element in a time step. Following this step, the newly convected element concentra- 
tions are dispersed laterally to each of the adjacent elements. Then the growth or the 
decay takes place. And all these three processes are simultaneous in nature. In each 
individual element it can be observed that the pollutant concentrations advected to 
that element go on dispersing and decaying exponentially in each time step. The 
computer code developed in this study for simulating the transport process is given 
in Appendix A. 



Chapter 3 

Solution Results and Discussion 


3.1 Introduction 

This chapter is devoted to the discussion of the solution results. The flow simulation 
model assuming steady state GVF in an open channel and the pollutant transport 
simulation model was applied to an open channel with immobile bed. The cross 
section of the channel and other details are obtained from the data related to an 
experimental channel as reported in Luk et al.(1990). The motivation for the same 
is that experimental observation of actual pollutant movement in a similar channel, 
however, considering 2 — D flow in a slightly curved channel has been reported earlier 
(Luk et al.,1990). Although the simulation results reported here are for a straight 
channel and considering only 1 — D channel flow and 2 — D pollutant transport, 
those experimental observations provide a basis for comparision of the simulation 
results. 


3.2 Data Used for Simulation 

An artificial straight channel with irregular cross section is considered for both flow 
and transport model simulation. The cross sections are similar to those presented 
in Luk et al.(1990) for a mildly curved channel. The cross section at upstream end 
{X = Om) is assumed identical to that at X = 4m, from the upstream end for the 
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simulation purpose, in this study. 

In order to relate depth of flow with top width and the flow cross section area, 
depths of the channel upto 6.25cm at X = Om andX = 4m; 8.75 cm at X = 8m; and 
9.75cm at X = 12m, respectively were considered. The irregular cross sections of the 
stream channel, that is, X = Om; X = 4m; X = 8m; X = 12m are shown in FigZ.l. 
The experimental channel as reported in Luk et al.(1990) followed a sine curve with 
a wavelength of 1.9m and amplitude of 0.29m. The channel was 0.3m wide, with 
irregular bottom features. The channel was located in a 3mX20m sand basin with 
a sloping bed and sheet metal side walls. A conatant flow was discharged into the 
channel to scour the bed features naturally. Sand was fed in at the upstream end 
at a rate as close as possible to the erosion rate. After equilibrium was established, 
the water was drained and the features were preserved by spraying a dilute resin 
solution on the bed. It was hardened with a dilute solution of formic acid resulting 
in a hard bottom. For the simulation in this study, the same cross sections were 
considered, with a straight longitudinal profile. The slope and discharge remain the 
same. 

In the flow simulation model the Manning’s coefficient of roughness, n is assumed 
to be 0.010 for an artificial channel. As the cross sections are given at an interval 
of 4m the same step size of 4m in the length direction is assumed. The paramaters 
used in the flow simulation model are presented in Table2>.l. 


Description 

Parameter 

Discharge in the channel, Q 

Accelaration due to gravity, g 

Channel bottom slope, 5o 

Kinetic-energy correction factor, a 

0.00396 m?/s 
9.81 m/s^ 

0.0035 

1.0 


Table 3.1 : Flow Parameter Vsilues 

In the transport simulation model the stream is divided into seven stream tubes. 
Each stream tube is 0.0428m wide with a number of varying length elements. 
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The length of each element is the product of local velocity and the time step. 
The flow simulation model is one dimensional resulting in no velocity variation along 
a given cross section. Therefore, the length of stream tubes remain same at a given 
cross section. The time step is assumed to be Isecond in the transport model. The 
input data for this model is presented in the Tabled. 2. 


Description 

Parameter 

Width of the stream channel 

0.3 m 

Flow rate, Q 

0.00396 mVs 

Averaged flow depth, h 

0.04 m 

Average flow velocity, u 

0.35 m/s 

Mean bed slope, Sq 

0.0035 

Number of stream tubes 

7 

Concentration in 3''‘^st.tube 

495.51 mg/L 

Concentration in T^st.tube 

1123.7 mg/L 

Concentration in 5‘^st.tube 

2317.33 mg/L 

Time step, At 

1 sec. 

Transverse mixing coefficient, E~ 

0.00021 m7s 

Rate of chemical reaction, a 

0.0 s-^ 


Table 3.2 : Pollutant Transport Parameter Values 

The transverse mixing coefficient, used in this model is assumed to be non 
zero since the transport model is 2 — D and there is transverse movement of the 
concentrations bacause of dispersion even though the average velocity across any 
cross section is uniform, as the flow simulation model is one — dimensional. The 
mixing coefficient can however, account for transverse velocity variation, unmodelled 
in the I — D simulation model. The flow simulation is assumed to be steady state, 
while the transport simulation model is transient. The value is taken identical 
to the value used for the 2 - D flow in the experimental channel, so that the mixing 
effects are atleast accounted for in the transport model. 
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A slug injection is an insertion of certain concentration instantaneously. A conti- 
nous injection can be considered as an aggregate of infinite number of slug injections 
at a very small time interval. In this study single or multiple slug injections are 
considered for simulation purpose. It is assumed that the slug is injected at time 
t = Isecond in three adjacent stream tubes on either side of the central stream 
tube, that is.S'’'^, and 5^^ stream tubes to simulate better mixing affects. This 
is because in the modelled one - dimensional flow, complete mixing cannot take 
place until the pollutant has travelled long distances downstream from the point of 
release. 


3,3 Simulation Results 

The computations are carried out with an initial maximum depth of 2.65cm in the 
upstream cross section {X = Om) measured from the bottom most point of the 
stream channel bed. The average velocity is higher at this cross section since the 
area of flow is lower for steady flow conditions to be maintained. In the experimental 
results (Luk et ah, 1990), the depth profile and channel cross section at A' = Om is 
not reported. Therefore, cross sections are assumed to be identical at X = Om and 
A” = 4m , in this study. The depth of 2.65cm at A'’ = Om is obtained by iterative 
trial and error solution of the flow simulation model so as to achieve depth profiles 
at X = 8m and X — 12m as obtained for the experimental channel. 

Since the model is one — dimensional a uniform average velocity is simulated at 
any particular cross section. The depth profiles are well matched with the experi- 
mental data at the downstream sections. The maximum depth of flow at X = 8m is 
observed as 6.788cm. This differs from the experimental value by less than Ipercent. 
A depth of flow at X — 12m is obtained as 7.760cm. This depth of flow differs by 
a still lesser value, that is Q.bpercent. These results clearly indicate that the simu- 
lation model results are well matched on the downstream sections. In this model a 
decreasing average velocity can be observed along the direction of flow due to the in- 
creasing cross sectional area along the direction of flow. These results are intutively 
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expected in natural channels also. 

In the pollutant transport model O.IL of a tracer solution with a concentration 
of o8000mg/L is injected in three adjacent stream tubes, that is, 3’"'^, 4^^ and 
stream tubes for better simulation. The results show that at X = 4m from the 
upstream side, complete sectional mixing is not achieved and hence the peaks in 
the concentration profiles are higher in certain stream tubes in which the slug is 
inserted at X = Om. In the remaining stream tubes the peaks are comparitively 
lower. Further downstream, peaks are less concentrated in those stream stream 
tubes in which the pollutant was released upstream. Also, the pollutant spread in 
the transverse direction due to dispersion effects. The concentration profiles at a 
section X = 12m show further dispersion of the pollutant. The time-concentration 
profiles obtained for this simulation model as well as for both the MABOCOST and 
experimental models are shown in Fig3.2. A usual comparison of the concentration 
profiles (break through curves) show that in the combined use of 1 — D flow and 
2 — D transport model, assuming a straight longitudinal profile, the peaks are higher 
compared to MABOCOST and experimental results. However, again due to reduced 
mixing on dispersion effects, the time base of the break through curves are shorter in 
these simulations. It is also observed that the peaks are achieved faster in the central 
stream tubes. In the of-centre stream tubes peaks are lower than the MABOCOST 
and experimental results, also due to the lack of transverse velocity in the flow 
modelled. 

The computational results for estimating total mass of the pollutant at a cross 
section X = 4m is presented in the TableZ.Z. The pollutant concentration in 
an element at any section is the product of discharge and area under the time- 
concentration curve for that element. The discharge per second passing through an 
element is equal to its volume for steady flow conditions (as At = Isec.). The 
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Fig 3.2 Comparison of Simulation Results with Experimental and original MABOCOST Results | 

(contd 


Time (s) 



Time (s) 



X=8.0 m 

Z=0.19m 

1- D Sim. - 

2- D Sim. - 

Exp. 


- 

4- 

1 

L_. 

,\a 

/' \ + ++ 

/' ff 

+4111 1 1 








Fig 3.2 Comparison of Simulation Results with Experimental and original MABOCOSTResult 










25 


St. tube 

number 

(1) 

Area under Time-Cone. 

curve (mg-sec./L) 

(2) 

Discharge 

(L/sec.) 

( 3 ) 

Mass 

(mg) (2 X 3) 

( 4 ) 

1 

769.125 

0.311 

239.19 

2 

1136.85 

0.529 

601.39 

3 

1727.56 

0.407 

703.11 

4 

2708.36 

0.504 

1365.01 

5 

2653.91 

0.712 

1889.58 

6 

1375.47 

0.863 

1187.03 

7 

1367.90 

0.458 

626.44 



Total mass 

6611.80 


Table 3.3 : Mass Balance Calculations 

results show that the total mass of simulated pollutant concentration at this section 
differs by approximately 13% when compared with the input mass at X = Om. 
Also, this difference at X = 8m is observed as 11.87%. This mass balance error 
can be considered to be within the acceptable range for numerical simulation. The 
concentration profiles for all the stream tubes at X = 4m are given in Fig 3.3 

Another illustrative example with varying decay factors 0.1 and 0.05 show de- 
creasing concentrations downstream of the channel when compared with those ob- 
tained without considering any decay. These profiles are shown in FigZA. A change 
in the growth or decay factor affects the advected values to the neighbouring ele- 
ments in the flow direction. For conservative pollutants this factor is kept zero since 
no growth or decay takes place. 

The transport model is also tested for a case of two slugs injected at 10 seconds 
interval with same concentration and input as reported earlier. These results show 
two peaks in each profile with the second peak slightly higher than the first one 
because of a small residual concentrations left over at various sections at the time 
of second slug injection. The concentration profiles are shown in the FigZ.o. 
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Fig 3.3 Concentration Profiles for all Stream Tubes at X = 4m 

It is apparent from the previous discussion tliat 1 — D flow and 2 — D transport 
simulation gives comparable results with respect to the actual experimental data 
and can provide acceptable intial results before a more sophisticated 2 — D flow 
modelling is attempted. This is a viable method for pollutant simulation in wide 
open channels especially during dry seasons, when flow is nearly steady GVF, but 
w'ater quality is a critical factor. 
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Chapter 4 


Summary and Conclusions 


4.1 Summary 

In this study (i) A numerical model is developed for simulating 1 — D flow under 
steady state conditions for channels with irregular cross sections (ii) A simulation 
model is solved for simulating the 2 — D transport of a pollutant in the channel. 
This model is applicable for both 1 — D and 2 — D flow systems using the stream 
tube concept. 

The flow simulation model is formulated using the fourth order Runge Kutta 
method. The flow model is solved for an open channel with specified irregular 
cross sections assuming immobile bed conditions. The transport simulation model 
utilizes the flow simulation solutions to simulate pollutant transport in the channel, 
assuming advective and dispersive transport processes. 

The numerical solution is based on the modification of solution procedure sug- 
gested earlier for Mixing Anlysis Based on the Concept of Stream Tubes (MABO- 
COST) in a 2 — D environment. The transport model simulates the spatial and 
temporal distribution of pollutant concentration for given pollution sources. The 
simulated results are compared with reported experimental results for a similar 
2 — D system as well as reported 2 ~ D simulation results using M.ABOCOST. 
These comparisons and mass balance analysis show that the combined I — D flow 
and 2 — D pollutant transport simulation results may be used as approximate solu- 
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tions for the 2 — D system for preliminary assessment. Both conservative pollutants 
and pollutants undergoing decay process are considered. 


4.2 Conclusions 

This study shows the potential applicability of a 1 — £> open channel flow simulation 
model and a 2 — jD pollutant transport simulation model for obtaining preliminary 
solution for an open channel system with 2 — D flow. The feasibility of solving 
these two models in conjunction for conservative and non conservative pollutants 
are demonstrated using illustrative examples. Mass balance considerations validate 
the solution results for the illustrative problem. 

The 1 — D flow and 2 — D transport model gives comparable results vis — a — vis 
experimentally observed data for a 2 — D flow system. This conclusions may not be 
valid when the transverse velocities are significant. 
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APPENDIX A 

* Program for 2-D Pollutant Transport Simulation Model 

* 

* DECLARATION 

double precision tcC50,50, 100) ,fc(50,50,50) 
double precision wc(20) ,xc(20) ,yc(20) ,zc(20) 
double precision wh(20) ,xh(20) ,yb(20) ,zh(20) 
double precision wu(20) ,xu(20) ,yu(20) ,zu(20) 
double precision zlC20 , 100) ,zr(20 , 100) 
double precision .1(20,100) ,u(20,100) 
double precision d(20, 100) ,dx(20, 100) 
double precision dnl(20,100) ,dnr(20,100) 
double precision tbeta(20, 100) 
double precision dz,z,l,r,ez,q,x 
double precision alpha, beta 
double precision slef,srt 
int eger t ime , dt , imax , j max 
integer is,je,lt 

* 


c ommon/f ir s t /h , u , d , dx 
common/second/dz,z,l,r,q 
common/ third/ time , dt , imax , jmax 
common/f ourth/zl , zr 
common/f if th/ dnl , dnr 



open (unit=9 , f ile="zero") 
open (unit=10 , file="f our" ) 
open(unit=ll,file="eiglit") 
open(unit=12,file="twelve'') 
open (unit=13,file=" input") 
open (unit=14 , f ile="conc") 
open (unit=]„5 , f ile="output ") 

* Enter th.e following parameters 

read ( 13 , *) imax-, t ime , dt , ez , q , dz , z , alpha , beta 
do 10 i=l,ifflax: 
read(9,*)wc(i) ,wh(i) ,wu(i) 
read(10,*)xc(i) ,xh(i) ,xu(i) 
read(ll,*)yc(i) ,yh(i) ,yu(i) 
read(12,*)zc(i) ,zh(i) ,zu(i) 

10 continue 

x=0.0 
ltn=l 

21 do 20 ib=l,imax 

if ( (x . ge . 0 . 0) . and . (x , It . 4 . 0) ) then 
tc (ib , Itn , 1) =wc (ib) + ( (x-0) /4 . 0) * (xc (ib) -wc (ib) ) 
h ( ib , Itn) =wh ( ib) + ( (x-0) /4 . 0) * (xh ( ib) -wh ( ib) ) 
u ( ib , Itn) =wu ( ib) + ( (x-0) /4.0)*(xu(ib) -wu ( ib) ) 
elseif ((x.ge.4.0) .and. (x. It . 8 . 0) )then 
t c (ib , Itn , 1) =xc ( ib) + ( (x-4) /4 . 0) * (yc ( ib) -xc ( ib) ) 
h ( ib , Itn) =xh (ib) + ( (x-4) /4 . 0) * (yh ( ib) -xh( ib) ) 
u ( ib , Itn) =xu ( ib) + ( (x-4) /4 . 0 ) * (yu (ib) -xu ( ib) ) 
elseif ( (x. ge . 8.0) . and. (x. le . 12 . 0) ) then 
tc(ib,ltn,l)=yc(ib)+((x-8)/4.0)*(zc(ib)-yc(ib)) 
h ( ib , Itn) =yh ( ib) + ( (x-8) /4 . 0) * (zh (ib) -yh ( ib) ) 
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dx ( ib , Itn) =u (ib , Itn) *dt 

d (ib , Itn) = ( (h (ib , Itn) **2) * (u (ib , Itn) **2) *ez) / (q**2) 

tbet a ( ib , Itn) =dx ( ib , Itn) *dz*li ( ib , Itn) 

if (x . le . 12 . 0) then 

jmax=ltn 

endif 

20 continue 

if (x . gt . 12 . 0) then 

goto 22 

else 

write(*,*)ltn,x 
x=x+dx(l,ltn) 
ltn=ltn+l 
goto 21 
endif 

22 write(*,*) ’At what location no. of st.tub,elment ’ 

read(=i' , *)mst ,mel 
read(14,*) (tc(i, j ,mt) ,i=3,5) 
it=l 
num=l 
goto 25 

15 do 30 ii=l,imax 

tc(ii ,num, it)=tc(ii ,num-l,it) 

30 continue 

25 do 40 is=l,iinax 

je=nuiii 

do 35 lt=it,time,dt 

call left(tc,is, je,lt ,slef) 

call right(tc,is, je,lt,srt) 

tc(is, je,lt+dt)=tc(is , je,lt)+(dt/theta(is, je))*(srt+slef) 
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tc (is , je j It+dt)— tc (is , je , It+dt) *exp(alplia*dt) 

35 continue 

40 continue 

it=it+dt 
nuin=num+l 

if (num.gt . jmax)goto 45 
goto 15 

45 do 50 kt=l ,time,dt 

write (15,100)kt,tc (mst ,mel , kt ) 

100 format(2x, 'Time' ,2x,i5,2x, 'Concentration' ,2x,f8.2) 
50 continue 

stop 
end 

subrout ine left (tc,is,je,lt,slef) 


* DECLARATION 

double precision tc(50,50,100) ,fc(50,50,50) 
double precision 21(20,100) ,zr (20, 100) 
double precision h(20,100) ,u(20,100) 
double precision d(20, 100) ,dx(20,100) 
double precision dele (20, 100) 
double precision dnl(20,100) ,dnr(20,100) 
double precision dz,z,l,r,q 
double precision slef,srt 
integer time,dt,imax, jmax,p 
integer is,je,lt 
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common/f irst/h , u , d , dx 
c ommon/ s e c ond/dz , z , 1 , r , q 
c ommon/ third/t ime , dt , imax , j max 
c ommon/f ourth/zl , zr 
c ommon/ f if th/ dnl , dnr 

* 


do 200 mi=l,imax 

do 220 mj=l,jmax 

if (mi . le . Dthen 

delc(mi,mj)=tc(mi,mj ,lt) 

dnl (mi , m j ) =h (mi , m j ) *u (mi , m j ) *dz/q 

else 

gl=0.0 

do 225 isk=l,mi 
gl=gl+h(isk,mj) 

225 continue 
gm=0 . 0 
mik=mi-l 
do 230 isl=l,mik 
gm=gm+h(isl,mj) 

230 continue 

dnl (mi , m j ) = (gm-gl) *u (mi , mj ) *dz/q 
delc(mi,mj)=(tc(mi-l,mj ,lt)-tc(mi,mj ,lt)) 
endif 

220 continue 
200 continue 
1=1 


m=l 
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do 250 p=l,m 
if (is .le. Dthen 
slef=0.0 
else 

avgl=(h(is, je)*d(is,je)+li(is-l, je)+d(is-l,je))/ 2.0 

prodl=l*dx(is-l , je) *avgl*dz*delc (is , je) 
f t=prodl/ (dnl (is , j e) **2) 

ftemp=((li(is-l, je)*dz*d(is-l, je))-(h(is, je)*dz*d(is, je))) 

ftemp=abs (ftemp) 

sd= (l*dx (is-1 , je) *f temp*delc (is , j e) ) / (dnl (is , je) **2) 

slef=ft+sd 

endif 

250 continue 
return 
end 

subroutine right (tc , is , j e , It , srt ) 


* 


* DECLARATION 


double precision tc(50,50,100) ,fc(50,50,50) 
double precision zl(20 , 100) ,zr(20 , 100) 
double precision h(20,100) ,u(20,100) 
double precision d(20,100) ,dx(20,100) 
double precision rdc(20,100) 
double precision dnl (20, 100) ,dnr (20, 100) 
double precision dz,z,l,r,q 
double precision slef,srt 
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integer time,d.t,imcLX, jmaxjp 
integer is,je,lt 


common/f irst/h. , u , d , dx 
conunon/second/dz , z , 1 , r , q 
c ommon/ third/t ime , dt , imax , j maix 
common/f ourth./ zl , zr 
common/fifth/ dnl , dnr 


do 300 il=l,imax 

do 320 jl=l,jmax 

if (il . ge . imax) then 

rdc(il, jl)=tc(il, jl,lt) 

dnr(il, jl)=h.(il, jl)*u(il, jl)*dz/q 

else 

rdc(il, jl)=(tc(il+l, jl,lt)-tc(il, jl,lt)) 
vb=0.0 

do 325 jsk=l,il 
vb=vb+h(jsk, jl) 

325 continue 
vc=0 .0 
ilm=il+l 

do 330 jsl=l,ilm 
vc=vc+b(jsl, jl) 
continue 

dnr(il, jl)=(vc-vb) *u(il, jl)*dz/q 
endif 


330 



300 continue 


r=l 

k=l 

do 350 p=l,k 

if ( i s . ge . imcLx) then 

srt=0 . 0 

else 

ravg=(h(is, je)*d(is, je)+h(is+l, je)*d(is+l, je))/2.0 

rmult=r*dx (is+1 , j e) *dz*rdc (is , j e) *ravg 
rtd=rmult/ (dnr(is , je) **2) 

serve=((h(is+l,je)*dz*d(is+l,je))-(h(is, je)*dz*d(is, je))) 
serve=abs (serve) 

f th= (r*dx (is+1 , j e) *serve*rdc (is , j e) ) / (dnr (is , j e) ** 2 ) 

srt=rtd+fth 

endif 

350 continue 
return 
end 

* * * 

TYPICAL INPUT TABLE : 

7 100 1 0.00021 0.00396 0.04285 0.0 0.0 

495.51 

1123.8 


2317.33 



